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INTRODUCTION 


Asymptotic analysis was the predominant method for the approximate 
analysis of shell structures before the advent of the computer. The Geckeler 
approximation, ref. (1) , for the solution of axisymmetric problems for shells 
of revoluton is based on asymptotic analysis, ref. (2). The purpose of this 
memorandum is to review asymptotic methods of analysis to identify features 
that can be adapted to numerical solutions for shell problems. 

Most of the literature on asymptotic solutions is based on linear shell 
theory, ref. (3). In linear theory, the coefficients of the differential 
equations of the theory can be written in closed form for a broad class of 
problems. The classical asymptotic analysis expands these coefficients in 
power series. The review of asymptotic methods here focuses on eliminating 
explicit expansions of the coefficients. Current applications require 
solutions of nonlinear problems in shell structures. The nonlinear 
differential equations of shell theory can be solved by Newton's method, 
which reduces a problem to a sequence of linear differential equations . The 
linear equations have coefficients that vary with the independent variables 
and with the current approximate problem solution. For viable numerical 
solutions, formal series expansions of the coefficients is not feasible. 

The development of asymptotic methods in shell analysis paralleled the 
analysis of the special functions of mathematical physics. Going back to the 
same starting point, this memorandum is concerned with the derivation of an 
algorithm for solving second order linear differential equations. The 
numerical algorithm is applied to several examples for comparison with 
analytical methods. Bessel's equation is examined as an example of an 
equation with an irregular singular point at infinity, ref. (4). The results 
of classical analysis, used in subroutines to compute Bessel functions of 
large arguments, are compared to the numerical algorithm. The Falkner-Skan 
equation from boundary layer theory is solved as a nonlinear example using 
Newton's method, ref. (5), and the numerical algorithm derived in this 
memorandum. 
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FIRST ORDER DIFFERENTIAL EQUATIONS 


The basic numerical algorithm for second order equations is based on the 
application of Newton's method to first order nonlinear equations. The 
first order equation is the simplest example for solving differential 
equations by Newton's method. Its solution is outlined in general form in 
this section before applying the algorithm to second order equations in the 
next section. 

Newton's method approximates the solution y=y(x) for the general 
equation, 

^ +F(x,y)=0 (1.1) 


by solving the following sequence of linear equations: 



+ P 


(m-1) 


(x) 



E ( m " 1 ) 


(x) 


m-1, 2 ,3,... (1.2a) 


(m) (m-1) f (m) 

y - y + 6 y 

_,(m-l) „ . (m-l) N 

P - F (x,y N ') 

* y 


(1.2b) 


(1.2c) 


(1. 2d) 


In equations (1.2), the integer m is an iteration counter. After m 

iterations, the function y^ is the approximation to y, the solution of the 

nonlinear problem, equation (1.1). The function £y^ m \ defined by equation 
(1.2b) as the correction to the solution after (m-1) iteration cycles, is 
computed in the mth iteration by solving the linear variational equation, 

equation (1.2a). The right-hand side of equation (1.2d) identifies E^ m ^(x) 
as the residual error when the (m-l)st approximation is substituted in the 
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nonlinear equation. Finally, the subscript comma followed by y on the right- 

hand side of equation (1.2c) denotes P^ m ^as the partial derivative of 
F(x,y) with respect to the dependent variable. The partial derivative 
appears because the linear variational equation results from expanding the 
nonlinear equation in a Taylor series about the (m-l)st approximate solution 
for y while holding x constant and from truncating the series after the 
linear term. 

* 

At each iteration step , the linear variational equation, equation 
(1.2a), is a first order equation whose general solution can be written 


6 y - C e 


-u(x) 


+ e 


-u(x 




u(s) 


E(s) ds 


(1.3) 


where 

u(x)- | P(t)dt (1.4) 

For specific problems, the integrals in equation (1.3) can be written in 
closed form or evaluated by approximate methods. One approximate method for 
integrals is asymptotic expansions that are closely connected with asymptotic 
expansions for differential equations. Rather than approximating specific 
integrals as they arise in Newton's method, the alternative adopted here is 
to replace the integrals with numerical quadrature formulas . 

QUADRATURE FORMULAS 

For numerical solutions by integration, accuracy is improved by avoiding 
numerical computation of the derivative appearing in the residual error, E. 
Combining equation (1.2a) and equation (1.2b) to eliminate Sy gives 

d y (m) j. „<m-l), x r (®) (m-1), (m-1). 

ax + p (*) [y - y ] - -F(x,y ') 


* 

The superscript (m) for the iteration will be dropped in equations 
that follow when no confusion arises. 
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or 


^ (m) + p(»-D 


(X) 


y (m)_ Q (m- 1) 


(x) 


(1.5) 


where 


^(m-1), . _(ra-l), . (m-1) . (m-1). 

Q (x)- P (x) y - F(x,y ) 


( 1 . 6 ) 


It is also useful in numerical solutions to have the capability to 
integrate both forward and backward. Formulas for integrating forward from 
x=b and for integrating backward from x=c are : 

u(x) - [ P(s) ds x > b (1.7a) 

J b 


u(x) - - f P(s) ds 
“ x 


X < c 


(1.7b) 


y(x) - y(b) e * U ^ X ^ + I 1 (b,x) 


x > b 

SM 


(1.8a) 


I^(b,x) - J e Q(t) dt x > b 


(1.8b) 


y(x) - y(c) e " u ^ x ^ - I 2 (x,c) 


X < c 


(1.8c) 


I 2 ( x 


,c) - |° e -l u <*>-u<t 


)] 


Q(t) dt x < c 


(1 • 8d) 


When P(x) is a smooth function, the integrals for computing u(x) , 
equations (1.7), can use Simpson's rule or other standard quadrature 
formulas. However, when P(x) is slowly varying but large in absolute value, 
asymptotic analysis of the convolution integrals 1^ and in equations 
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(1.8), ref. (4), suggests using the exponential function, exp { - [u(x) -u( t) ] } , 
as a weighting function in the quadrature formula. The exact form of the 
exponential functions is not known in advance, but recurrence relations for 
the integrals allow use of approximate weighting functions. 

RECURRENCE RELATIONS 

The convolution integrals for the particular solutions obey the 
following recurrence relations between integrals for two consecutive values 
of x, say x-x^ and x-x^ + ^: 

I l^ b,X i+l ) “ exp( - [u(x i+1 ) -uCx^ ] ) + I 1 ( x i > x i+1 ) (l-9a) 

I 2 (x i( c) - I 2 ( x i + i' c ) exp{-[u(x i )-u(x i+1 )]) + I 2 ( x i> x i+i) (l-9b) 

For the numerical quadrature formulas , the interval of integration is 
divided into subintervals determined by a set of stations at x-x^. From the 

definition of u(x) , equation (1.4), and the mean value theorem, 

u(x) - u(t) - P(£) (x-t) t < £ < x (1.10) 

The mean value theorem suggests that a weighting function of an exponential 
term with a linear argument can be used for each sub interval. Let 

g(t) - u(x) - u(t) - a(x-t) a-P(x) (1.11) 

Then, changing the dummy variable reduces the integrals in each subinterval 
to the standard forms , 


h^i-v V- L 1 


e aS f^(s) ds 


v. - x. - x. 1 (1.12a) 

l l l-l 


V X i ,X i+l )= J Q 1+1 eHS f 2 (s) ds a = p ( x i) 


(1.12b) 
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where 


f^s) - exp [-g(x i *s) Q(x^-s) ] (1.12c) 

f 2 ( s ) - exp [ - g (x^+s ) Q(x^+s) ] (1 . 12d) 

The final steps in deriving quadrature formulas for the integrals in 
equations (1.12) are to expand the functions f(s) with polynomial 
interpolation formulas from tabulated values of P, Q, and u and to integrate 
term by term. A useful formula for summing the integrals is 




ds 


n -as 
s e 

a 



(1.13) 
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SECOND ORDER LINEAR EQUATIONS 


The results derived in the preceeding section for first order nonlinear 
differential equations are readily extended to nonlinear second order 
equations, ref. (5). The linear variational equations of Newton's method at 
any iteration cycle will have the general form 

y" + p(x) y' + q(x) y - - e(x) (2.1) 

where primes denote differentiation with respect to x. In Newton's method, 
the dependent variable y is a correction to a current approximation. The 
concern here is not with the derivation of equation (2.1) for a specific 
problem, but with its numerical solution. The numerical solution does not 
operate on equation (2.1) directly, but only a related nonlinear form, 
suggested by asymptotic theory. 

Asymptotic analysis starts by assuming the complementary solution of the 
second order linear equation in the form 

y(x) - exp [J <f>( t) dt] (2.2) 


For a solution of the complementary equation, <£(x), the derivative of the 
integral at x, must satisfy the Riccati equation, ref. (6), 

<f> • +^ 2 + p^ + q- 0 (2.3) 

The assumed solution for y(x) , equation (2.2), transforms the linear 
second order equation into a first order nonlinear equation for <£(x). The 
first order equation will be solved here by Newton's method as outlined in 
the section on first order equations. The linear variational equation, 
equation (1.2a), in an iteration for <f> is 


(5^') (m) +[2 p] - - E[^(x)] (m ‘ 1) 


(2.4) 
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and equation (1.5) for the Riccati equation is 


(^') (m) + [2 * (m_1) + p] / m) - [(* 2 ) (m_1> - q] (2.5) 

The iteration by Newton's method requires a zeroth approximation to start the 
iteration. For numerical solutions, a boundary condition for <j> is also 
required. Insight into stable numerical solutions of the Riccati equation is 
found in the theory for second order linear equations . 

LINEAR INDEPENDENCE OF SOLUTIONS 

A fundamental set of solutions for a linear second order differential 
equation is defined as a pair of linearly independent solutions. Therefore, 
two solutions, <f > ^ and 4 > must be found for the Riccati equation so that a 

fundamental set of complementary solutions for y can be written in the form 

y “ C 1 y x (x) + c 2 y 2 (x) (2.6) 

with 

y £ (x) - exp [u ± (x) ] (2.7) 

u i “ | <^ i (t) dt i- 1,2 

To be linearly independent, the Wronskian, W, of the solution must not 
vanish. The Wronskian of a fundamental set of solutions is 

W - ( y 1 y 2 - y 2 y{) (2.8) 

For the fundamental set, equations (2.7), defined by solutions of the Riccati 
equation, the Wronskian is 


P x 

W(x) - [* 2 (x) - ^(x)] exp { [<f> 2 (t) + ^(t)] dt } 


(2.9) 
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Since the exponential function never vanishes, the solutions are a 
fundamental set as long as <f>^ and <j > 2 are distinct solutions of the Riccati 

equation. Therefore, the choice of zeroth approximation and boundary 
condition in solving the Riccati equation, equation (2.3), by Newton's method 
must be influenced by the requirement of finding two solutions that do not 
intersect in the interval of integration on x, the independent variable. 

COMPUTATION OF WRONSKIAN 

Solving the linear variational equations by computing convolution 
integrals numerically has the potential for accurate and numerically stable 
solutions. The Wronskian is also a measure of the accuracy of the two 
solutions, <f > ^ and <f>^. The difference in the two solutions can be checked 

numerically by a direct integration rather than a convolution integral . The 
difference can be derived by substituting <t> ^ and <f>^ in the Riccati equation 

separately and subtracting the two equations . The resulting differential 
equation for the difference (<f > 2 - <j>^) is 

(<f>2 * h y + (p + *2 + *1> ( V V " 0 (2.10) 

The solution of equation (2.10) can be written as 

yi 2 (x)-^ 1 (x) - [^(b)-^ 1 (b)] exp {- Jjp(t)+^ 2 (t)+^ 1 (t)] dt) (2.11) 


where x-b is any point in the interval of integration. The equation for the 
difference in the two solutions also follows from substituting Abel's 
identity for the Wronskian, 


W(x) 


C exp( 


•r 


P dt) 


( 2 . 12 ) 


into equation (2.9). 
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ZEROTH APPROXIMATION AND BOUNDARY CONDITIONS FOR <f> 


The theory for the second order linear differential equation requires 
two distinct solutions for the Riccati equation. The solution for the 
Riccati equation by Newton's method is influenced by the zeroth approximation 
and by the constant of integration for the linear variational equation in the 
iteration. When the coefficients p and q in the second order equation are 
slowly varying functions of x, then classical asymptotic theory suggest two 
distinct zeroth approximations. First, the theory changes the dependent 
variable in the Riccati equation to a, where 

a = (j> + (p/2) (2.13) 

The change of variable inserted in the Riccati equation, equation (2.3), 
transforms it to 

a' + a 2 + 5 - 0 (2.14) 


where 


q - q - (P 2 /M -(p'/2) (2.15) 

Then, two zeroth approximations to solve the transformed Riccati equation by 
Newton's method are 


a 


( 0 ) 

1 



(2.16a) 


a 


( 0 ) 

2 




(2.16b) 


When q - q(x) is a real function, the two zeroth approximations for a are 
both real or both pure imaginary. For the case with complex a, only one 
iteration sequence is required starting with one of the two zeroth 
approximations. If a is complex, then by definition, equation (2.13), <f> is 
also a complex solution of the Riccati equation, equation (2.3). It follows 
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by direct substitution that if <j > ^ is a complex solution, then ^2 ” ^1* t * ie 
complex conjugate of <j > is another solution. 


The second case where q<0 and the two zeroth approximations for a are 
real requires two iteration sequences in Newton's method starting each 
sequence with one of the zeroth approximations. If the two converged 
solutions, and , do not vanish in the interval of integration and retain 

the algebraic signs of the zeroth approximations, their difference, a < 2 ^, 

will not vanish. By definition, equation (2.13), the difference (<f>^ - * 

“ c*^) will not vanish either. The pair of solutions of the Riccati 

equation, <f>^ and <f > when integrated determine two exponential functions 
that are two linearly independent solutions for y in equation (2.6). 

In either case, when the imaginary part of a complex <j> does not vanish 

or when two real a of opposite sign are found, the Wronskian will not vanish 

because of equation (2.9) and the exponential solutions for y will be a 

fundamental pair of solutions for the linear second order differential 

equation, equation (2.1). 

BESSEL'S EQUATION OF ORDER n 

An example problem will clarify the procedure for solving the Riccati 
equation by Newton's method for two distinct solutions. Consider Bessel's 
equation of order n for positive values of the real variable x, 

y" +(l/x) y' + (1-nV x^) y - 0 0< x < ® (2.17) 

For this problem p = (1/x) and 

q - [l-(n 2 - l/4)/x 2 ] (2.18) 

For this example, the algebraic sign of q depends on the value of x 
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q > 0 


(2.19) 


x > £ = (n 2 - 1/4) 1/2 


The point x - £ is called a turning point or a transition point, because of 
the change in qualitative behavior in y as the exponential solutions , 
equation (2.6), pass from real to imaginary as x passes through the point. 
The classical literature seeks a uniform asymptotic expansion valid on both 
sides of the turning point. For numerical solutions, the additional analysis 
is not necessary since linearly independent solutions on either side of the 
turning point can be matched by a direct computation to determine constants 
of integration. Therefore, solutions to the left and right of the turning 
point for Bessel's equation are computed separately here. 


Consider first, solutions of the Riccati equation for Bessel's equation 
in the interval b < x < £ where b > 0 to exclude the regular singular point 

at x — 0 . The zeroth approximation 


a 


( 0 ) 

1 


2 2 

(r/x - 1) 


1/2 


for n-1 is shown in figure 1 along with the approximation after one iteration 

cycle of Newton's method, and the converged solution, . The first 

approximation lies above the zeroth approximation and is positive at the 
turning point of the zeroth approximation. 

The first approximation for a remaining positive when the zeroth 
approximation is positive is a general result. Because the zeroth 

approximation satisfies equation (2.16a), the residual error in the linear 
variational equation for the first iteration is minus the derivative of the 
zeroth approximation 

8al^ +2a^ (2.20) 

where superscripts denoting the first iteration have been dropped. 
Integrating the variational equation between limits gives 
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<-ap dt 


( 2 . 21 ) 


Sa 1 (x) - [ e 
J b 


-(u(x)-u(t)] 


Since 


u(x)-u(t) 



2a^(s) ds > 0 


( 2 . 22 ) 


the exponential function in the integral in equation (2.21) is less than 
unity, 


0 < e -[u(x)-u(t)] 


< 1 


(2.23) 


For cases where < 0 in the interval of integration, as is true for the 

zeroth approximation in figure 1, the properties of the integral in equation 
(2.21) shows that the first correction to is always positive, 


0 < £a^(x) 



dt - a^(b)-a^(x) > 0 


(2.24) 


or writing the result of the first iteration with superscripts, 


6a™ + a<°> > a<°>> 0 


(2.25) 


When the slope of the zeroth approximation is positive in the interval of 
integration rather than negative, the particular solution for the first 
correction will give the same result as equation (2.24) by interchanging the 
end points to integrate from x to b. 

It also follows by a similar proof that the second correction in 
Newton's method is negative. The converged solution for a ^ in figure 1 lies 

between the first approximation and the zeroth approximation, but the proof 
that this is true is problem dependent and not examined here. Also, the 
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numerical solution does not use the linear variational equation directly to 
solve for the correction Sol because the derivative of the zeroth 
approximation does not exist at the turning point £. The form corresponding 
to the general equation (1.5) is used for the numerical solution for a. 

Once the positive solution, a^, of the transformed Riccati equation, 

equation (2.14), for Bessel's equation is computed, a second iteration 
sequence of Newton's method is required to compute a second solution, a^- 

The result of this iteration, starting with a negative zeroth approximation, 
is shown in figure 2. The integration direction in the particular solutions 
for this case is backward from the turning point, x=£. The backward 
integration makes the first correction lie between the x-axis and the zeroth 
approximation. The converged solution for a 2 is less than or equal to zero, 

is always greater than zero. In the interval 0 < x < f, Newton's method 

gives two real a of opposite sign and two linearly independent solutions of 
Bessel's equation follow from computing two solutions of the original Riccati 
equation, 

<f> ± - ( -l/2x) + i-1,2 (2.26) 

The solution of Bessel's equation is then written in the form 

y - exp[u 1 (x)] + C 2 exp[u 2 (x)] 0 < x < £ (2.27) 

where 

u^x) - - f ^(t) dt < 0 
J x * 

u 2 (x) - | <^ 2 (t) dt < 0 

The directions of integration to compute the integrals u^(x) are chosen to 

normalize the exponential solutions of Bessel's equation to less than or 
equal to unity. 
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COMPLEX EXPONENTIAL SOLUTIONS OF BESSEL'S EQUATION 


For solutions of Bessel's equation to the right of the turning point, 
x > £, the two zeroth approximations for a in equation (2.16) are pure 

imaginary complex conjugates. In this case, only one solution of the 
transformed Riccati equation is required for a as long as the imaginary part 
of the complex solution does not vanish. Since complex algebra is required 
in the numerical solution, the amount of computation is approximately the 
same for one complex solution for a as for two real solutions . The 
quadrature formulas derived for first order linear real equations are 
identical for complex equations so that computer subroutines in real 
variables are readily converted to complex subroutines . 

The derivative of the zeroth approximation for a approaches zero as x 
goes to infinity so that the residual error in the linear variational 
equation, equation (2.20), is small for large values of x. Setting the 
correction equal to zero at a finite x and integrating backward in the 
quadrature formulas gives rapid convergence for Newton's method to compute a 
complex solution a of the transformed Riccati equation. Real and imaginary 
components of for the case of order n-1 are shown in figure 3 . 


With — a^, a complex solution for Bessel's equation follows from 

equation (2.26) and equation (2.27), or the final result for y can be written 

in real form since C2~ C^. The numerical results for y from Newton's method 

compare favorably with tabulated values for Bessel functions of the first and 

second kinds , J (x) and Y (x) . The quadrature formulas for the interval of 
n n 

integration in figure 3 used 41 stations. The numerical solution for J^(x) 

and Y^(x) is accurate to six digits. The solutions for the Bessel functions 

are less accurate to the left of the turning point reflecting the rapid 

variation in a caused by the regular singular point at x=*0. 
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IMPROVED ANALYTICAL RESULTS FOR BESSEL FUNCTIONS 


The numerical results for the complex solution of the tranformed Riccati 
equation, equation (2.14), for Bessel's equation suggest an analytical 
improvement in the zeroth approximation for a in the general case. The real 
component of a in figure 3 is small compared to the imaginary component. A 
look at the first iteration in Newton's method suggests that a small real 

part for a complex a can be expected when q is positive, large, and slowly 
varying. 

Dropping the superscript zero in equation (2.16a), 


a 


(q) 1/2 


Also, differentiating 


2 

a 


-q gives the expression for the derivative, 


a' - -q' /2a (2.28) 

The derivative is the residual error appearing in the first iteration of 
Newton's method for the correction 5a, 


5a' + 2a 5a - -a' (2.29) 

When a' is small in absolute value compared to a and slowly varying, an 
approximate particular solution of equation (2.29) is obtained by neglecting 
the derivative on the left-hand side and substituting for a' from equation 
(2.28), 

5a - -q'/4a 2 -q'/4q (2.30) 

Therefore, adding the approximate correction to the zeroth approximation 
gives a first approximation, 

a (1) = (q ' / 4q) + i q 1/2 (2.31) 
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The corresponding approximation to the solution of the Riccati equation 
is 


<f> — -p/2 + a 


( 1 ) 


The approximate complex solution for y is 


y = (q) " 1//4 exp [ -J (p/2) dx + i J q ^dx] 


(2.32) 


For Bessel's equation, the approximate solution for y can be written in 
closed form. By selecting a complex constant of integration that makes the 
approximate complex solution agree with the Hankel function of the first kind 
as x goes to infinity, the real and imaginary parts of y approximate the real 
Bessel functions of the first and second kinds: 


J (x)= (2/ttX) 1/2 cos [X + £ tan _1 (£/X) - nn/2 - n/U] (2.33a) 

Y n (x)- (2/irX) 1/2 sin [X + £ tan _1 (£/X) - ror/2 - ff/4] (2.33b) 


where 


2 2 2 2 2 

X - x - ^ and f - n - 1/4 

The corresponding result from the pure imaginary zeroth approximation, 
ct-iX/x, is 


J (x)= (2/irx) 1/2 cos [X + | tan‘ 1 (C/X) - nir/2 - w/4] (2.34a) 

Y n (x)= (2/irx) 1/2 sin [X + £ tan _1 (£/X) - m/2 - rc/W] (2.34b) 
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The approximation from Newton's method for the two solutions of Bessel's 
equation can be compared to the leading term from the classical asymptotic 
expansion, 

J^Cx)- (2 /ttx) 1 / 2 cos (x - xm/2 - 7r/4) (2.35a) 

Y^(x)- (2/ttx) 1//2 sin (x - n?r/2 - tt/4) (2.35b) 

or the first approximation from perturbation theory, ref. (7), 

Jn( x )-“ (2AX) 1 / 2 ! cos (x - n7r/2 - tt/4) 

- (£ 2 /2x) sin (x - n7r/2 - tt/4) ] (2.36a) 

Y n (x)- (2/ttx) 1/2 [ sin (x - mr/2 - tt/4) 

+ (£ 2 /2x) cos (x - mr/2 - tt/4) ] (2.36b) 

Numerical results from the various approximations for Jg(x) are compared to 

exact values in Table 1. The approximation in equation (2.33a) derived from 
the approximate first iteration for a, equation (2.31), in Newton's method is 
the most accurate of the approximations in the interval £ < x < ® except 

near x-£. The approximation breaks down at the turning point £ since the 
derivative on the left-hand side of the linear variational equation, equation 
(2.29), cannot be neglected near the turning point. 

The general case will give results similar to those for Bessel's 

equation. Whenever q is large, positive and varies slowly with x, equation 
(2.32) is a good approximation for a complex solution of the second order 
linear equation (2.1) 

INTEGRAL REPRESENTATION OF THE PARTICULAR SOLUTION 

The general solution of a second-order linear differential equation 
requires a particular solution in addition to two solutions of the 
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homogeneous equation. Variation of parameters is a general method for 
determining particular solutions by integration. The integrals that appear 
in the particular solution are of the same form as those for the solution of 
first order linear equations when the solutions, equations (2.7), of the 
homogeneous equations have the exponential form derived from the solutions of 
the Riccati equation. The general expression for the particular solution of 
the inhomogeneous equation, equation (2.1), is 


y p (x) 


y i (x ) J [y 2 (t) 


e(t)/ W] dt -y 2 (x)J* [y 1 (t)e(t)/W] dt 


(2.37) 


where W is the Wronskian of the two solutions, y^(t) and y 2 (t). When y^ and 
y 2 are of the exponential form in equation (2.7) and the Wronskian is also 
exponential, equation (2.9), the particular solution is 


y p (x) “ y pl (x) + y p2 (x) 


(2.38) 


where 

y p l< x ) “ J exp [u x (x) - u L (t)] (e(t)/[^ 2 (t) - -^(t)]} dt 

y p2 (x) “ ’ J 6Xp f u 2 (x) ' (e(t)/[^ 2 (t) - <^ 1 (t)]) dt 

The integrals in the particular solution are of the same type as the 
convolution integrals, equations (1.8), representing the solution of first 
order linear differential equations. The numerical quadrature formulas 
derived for an exponential weighting function, equations (1.12), apply to the 
particular solution with the <f > functions in the role of the weighting 

function, P(x). 

The properties of the method of variation of parameters make the 
derivative of the particular solution simple to compute from the expression, 

y p (x) = ^l (x) y pl (x) + ^ 2 (x) y p2 (x) (2.39) 
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ASYMPTOTIC INTEGRATION APPLIED TO THE FALKNER-SKAN EQUATION 


The numerical analysis of the previous 
Falkner-Skan equation. 

section is applied here 

f' " + f f" + £ [1 - (f ' ) 2 ] - 0 

(3.1) 

with boundary conditions , 


f(0) - 0 

(3.2a) 

f'(0) - 0 

(3.2b) 

i-b 

/--N 

8 

V-/ 

1 

(3.2c) 


The Falkner-Skan equation is a third-order nonlinear ordinary differential 
equation derived from Prandtl's partial differential equations for fluid flow 
in laminar boundary layers, ref. (8). In order to test the method of the 
previous section, the third- order system is solved here by a modified form of 
Newton's method. The modification reduces each iteration cycle to solving a 
linear second order equation and integrating a simple first-order equation. 

The modified form of Newton's method converges rapidly for positive 
values of the parameter and more slowly for negative values. The 
convergence is somewhat surprising since recent literature conveys the 
impression that Newton's method is unstable for the Falkner-Skan equation, 
ref. (9) . The analysis in this section shows that the instability is in 
shooting methods, ref. (10) and ref. (11), that use forward integration. The 
numerical instability is not in Newton's method per se. 

The convergence properties arise naturally as part of the iteration 
cycle. Two zeroth approximations that start the iteration and the resulting 
first approximations provide insight into the qualitative behavior of the 
final solutions. The qualitative behavior has been studied by Hartree, ref. 
(12) and by Smith, ref. (13), but they overlooked the approach of combining 
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asymptotic analysis and numerical analysis through Newton's method that is 
used here. 

THE TRANSFORMED NONLINEAR PROBLEM 

The zeroth approximations suggest a change in the dependent variable. 
The derivation of the iteration cyle begins with rewriting the problem in 
terms of w where 

f*x+w (3.3) 

The transformed nonlinear problem in w is 


w' ' ' + (x+w) w' ' -£(2+w') w'~ 0 


(3.4) 

$ 

/-N 

0 

s-/ 

1 

o 


(3.5a) 

✓—N 

0 

1 


(3.5b) 

w' (<*>) — 0 


(3.5c) 

MODIFIED NEWTON'S METHOD 



The linear form of Newton's method 
operator notation is 

for 

the mth iteration cycle 

L(«w<“>) - -e(w^ m *^ ) and „< m) - w (n - l> + 

5» <m) 

(3.6a) 

- L<w <m - 1) )-a(w <m - 1) ) 


(3.6b) 


The modified form of Newton's method does not update the coefficients in the 
linear operator L( ) with the current approximation for w but retains an 
earlier approximation w=w 0 . 
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Omitting details of the derivation, the modified form of Newton's method 
for solving equation (3.4) for w used here is 

y''+py'+qy * -e (3.7a) 


where 


y— w' 

(3.7b) 

p-(x+w 0 ) 

(3.7c) 

q— 20(l+w o ' ) 

(3 . 7d) 


e- (w-w 0 )w' r - /3(w'-2 w^)w' (3.7e) 

The system of equations (3.7) is solved iteratively by substituting the 
(m-l)st approximation for w into the expression for the residual e, equation 
(3.7e),and solving the linear system, equation (3.7a) and equation (3.7.b), 
for the mth approximation for w that satisfies the boundary conditions, 
equations (3.5). The iteration procedure is started by letting the zeroth 
approximation for w be equal to w 0 . The term e is not the usual residual 

error in Newton's method because of the L(w) term in its definition, which 
corresponds to the right-hand side of equation (3.6). However, when L(w) = 
0, then e is the usual residual error and will be referred to as the residual 
error below. 

ZEROTH APPROXIMATION AND FIRST ITERATION 

If w 0 - 0, the first iteration cycle solves the linear equations 


y' '+xy ' -2/3y = 0 


(3.8a) 


w' = y 


(3.8b) 
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The solution w gives a first approximation for the solution of the nonlinear 
problem in w, equations (3.3) and equations (3.4), and also a good 
approximation for f=x+w with the error in f equal to the error in w. 
Assuming the solution of the form, 

y - exp [ <f> dt (3.9) 


w' 


leads to the Riccati equation for equation (3.8a) 


<fi'+ f + x 4> - 2 p 


(3.10) 


As x goes to infinity, two distinct solutions of the Riccati equation 
can be identified, <f>^ - -x and ^ ~ 2/3/x • The solution for equation (3.8a) 

can be written as 


pX pX 

w' - c 1 exp [J ^dt] + c 2 exp [J t^dt] (3.11) 

where c^ and c 2 are constants of integration. For large values of x, an 
asymptotic expression for the solution is 

w' ~ c^exp (-x 2 / 2) + c 2 x 2 ^ 3 (3.12) 

When /9 is positive, the constant of integration c 2 must be zero in order 

to satisfy the boundary condition on w' at infinity, equation (3.5c). When /3 
is negative, w' vanishes at infinity for any real, finite choice of the 
constants of integration. Hartree, ref. (12), derived the equivalent of 

equations (3.8) and argued on physical grounds that the constant c 2 should 

also be set to zero for negative values of f3 . When l<2/3<0, it seems that a 
stronger argument could be made by noting that the integral of w' does not 
have a finite limit at infinity. The nondimens ional stream function f=x+w 
appears in the dimensional expression for the transverse component of the 
flow velocity, ref. (13), and the physics of the problem suggest requiring 
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that w be finite at infinity. In any case, the constant c 2 is suppressed in 
the numerical results reported here. 

Smith, ref. (13), observed that, when 2£ is zero or a positive integer, 
an exact solution of equations (3.8) can be written in terms of integrals of 
the error function. These integrals are computed for large values of the 
argument by asymptotic expansions. However, Smith did not pursue asymptotic 
expansions as a method of computation. Instead he used a power series 
solution for solving the nonlinear problem starting the solution at x=0 and 
using trial -and- error to establish the value of f ' ' (0) necessary to start the 
forward integration. 

The form of the asymptotic solutions for the first correction suggests 
the opposite approach of integrating backward from infinity. The Riccati 
equation can be solved numerically for two distinct solutions, <f>^ ~ -x and 

^2 ~ 2/9/x by integrating backward. For each iteration cycle, the numerical 

solution of equation (3.7a) finds the solution by the method of the previous 
section with the solutions of the associated Riccati equation computed by 
backward integration. This procedure identifies a constant of integration 
corresponding to that is expliticly set to zero for every iteration cycle. 

Before examining the numerical solution for the first approximation and for 
succeeding iterations for f, another choice for w 0 , the zeroth approximation 

for w, is introduced here. The second choice is w 0 - -c, where c is a 

constant. For this choice, the first iteration solves 

y' ' + (x-c) y' -2/9 y - 0 (3.13a) 

w' - y (3.13b) 

instead of equations (3.8), which becomes the special case of c—0. The 
choice of c is arbitrary. It is an input parameter that can be used to make 
the solution of equations (3.13) a better appoximation to the solution of the 
nonlinear problem in w, equation (3.4) and equations (3.5). The residual e 
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computed by substituting the solution of equations (3.13) into equation 
(3.7e) is a true residual error since equations (3.13) are homogeneous. 

e - (w+c) w' ' - f) (w') 2 (3.14) 

For positive values of /9, there is a real constant c that speeds the rate of 
decay of e as x goes to infinity. This choice is 


c - - [ w' (t) dt (3.15) 

J 0 

The reason for the improved rate of decay is clearer for the special 
case of £-n/2 where n is zero or a positive integer. For this case, the 
solution of equations (3.14) can be written in terms of functions with known 
asymptotic expansions. 

CHOICE OF c WHEN 2/9 IS AN INTEGER 

The solution of equations (3.13) for the first approximation for w can 
be written in terms of integrals of the error function when 

2/9 - n n- 0 , 1 , 2 , . . . 

For brevity, the usual notation, ref. (14), is abbreviated here. Let 

z - d (x-c) d 2 - 1/2 (3.16) 


,n 

l 


(z) 



i/n-1) 


(t) dt 


n- 0 , 1 , 2 , 3 , . . . 


(3.17) 


where 


-1 i /2 2 

i (z) - (2 /tt ' ) exp (-z ) 


and 
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o i /2 r 2 

i (z) = ( 2/7T ' ) exp (-t ) dt 

z 

The solution of equations (3.13) that satisfies the boundary conditions, 
equations (3.5), when 2fi-n is 


w' - -i n (z)/i n (-dc) (3.18a) 

w- [ i n+1 ( 2 ) - i n+1 (-dc) ]/[di n (-dc) ] (3.18b) 

For this case the constant c in equation (3.15) is defined as 

c -[ i n+1 (-dc) ] / [ d i n ( -dc) ] d 2 -l/2 (3.19) 

making 


w+c- i n+ ^(z)/i n ( -dc) - A exp (-z^)/ z n+ ^ (3.20) 

where A is a constant. The asymptotic behavior of w+c shows that the order 
of the residual error, e, equation (3.14), is quadratic in w' . 

The actual computation of e requires a numerical value for the constant 
c. The equation that defines the constant c for 2/3-n is transcendental, but 
an exact solution is not necessary since the first approximation for w is 
corrected later in the iteration. Graphical solutions, equating both sides 
of equation (3.19), are shown in figure 4 for different values of n. 
Interpolating the tabulated values used in plotting the curves in the figure 
gives the values for c for integer values of n that are listed in Table 2. A 
cross -plot for interpolating c when 20 is not an integer is shown in 
figure 5. 

Given c for each n, the residual error e can be computed using w and its 
derivatives from equations (3.18). The residual error is plotted as a 
function of x for integer values of n=2/3 in figure 6. There is no finite 
solution for c when n=0 in equation (3.19) and the residual error for n=0 is 
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from the error function with c-0. Even though the term in /3 is zero, the 
error for this case is larger for values of x away from the origin than the 
error for positive values of n associated with a finite value of c. 


NUMERICAL RESULTS FROM THE MODIFIED NEWTON'S METHOD 

The numerical solution of the Riccati equations for equations (3.13) has 
been programmed on a microcomputer for any input value of /3 and c. The 
solutions for fi - n/2 in terms of repeated integrals of the error function 
provide a check of the accuracy of the numerical solution. The next step in 
the general case, after solving equations (3.13), is to compute the residual 
e from equation (3.14) and start the iterative solution of equations (3.7). 
The inhomogeneous second-order differential equation, equation (3.7a), 
provides a test case for computing particular solutions from the integrals in 
equation (2.38). 

The first approximation for w from equations (3.13) is a special case of 
the first iteration of the general form of the modified Newton iteration, 
equations (3.7) where w 0 - -c and e-0. A second program was also written 

that updates w 0 with the current approximation for w from the first program 

that uses w 0 - -c. 

Numerical results from the two programs are summarized in Table 3 . All 
the numerical results were computed using 41 equally spaced stations for x. 
The numerical solutions for the special cases of repeated integrals of the 
error function are accurate to five or six digits. Based on published values 
for f''(0), the particular solutions computed from the integrals of equations 
(2.38) retain the same accuracy. 

The modified Newton's method gives good convergence for positive values 
of /3. The first program with w 0 - -c in the linear operator gives fair 

results for f ' ' (0) after a few iterations. The behavior at infinity is 
numerically stable so that the numerical solution follows the form of 
equation (3.12) plus a particular solution at each iteration step. The 
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constant c 2 **® for a finite solution at infinity and the constant c^=-l to 
satisfy the condition on w' at x-0. 

If the iteration of equations (3.7) is carried out using shooting 
methods with forward integration, suppressing the solution multiplied by the 
constant c ^ by adjusting f''(0) requires a good approximation for f ' ' (0) . 

The asymptotic analysis shows that a more stable shooting method would be to 
integrate backward with w(«>) as the initial condition to be determined. 

The modified Newton's method has poor convergence for negative values of 
0 ; the linear form of Newton's method also has slow convergence for this 
case, ref. (15). The slow convergence is a consequence of the lack of 
uniqueness of solutions of the nonlinear problem when 0 negative. 
Stewartson, ref. (16), showed the existence of solutions with negative values 
of f''(0) for certain values of 0 . Since the nonlinear problem is 
autonomous, it should be possible to alter the iteration along the lines of 
ref. (17) while retaining the numerical asymptotic integration to solve the 
linear variational equations. That approach requires more analysis and 
additional programming and is not explored further here. 
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CONCLUDING REMARKS 


An investigation was conducted to study numerical analysis of asymptotic 
expansions that arise in shell theory. The study showed that it is feasible 
to achieve the linear independence of the asymptotic theory for general cases 
while avoiding algebraic manipulations for each special case. Second-order 
nonlinear ordinary equations were examined in detail in this study. The 
results have desirable mathematical properties that can be extended to 
higher -order equations that appear in shell theory. 


The analysis of the Riccati equation derived from Bessel's equation 
produced a closed- form exponential approximation for Bessel functions of 
large argument that is more accurate than the leading terms of the usual 
expansion in negative powers of x. The result is general and applies to 
other problems whose solutions are nearly periodic over some range of the 
independent variable. 

Numerical quadrature formulas were derived to compute the convolution 
integrals that are part of the asymptotic analysis. The formulas use an 
exponential weighting function to produce a better fit of the integrand than 
a standard polynomial fit. 

The Falkner-Skan equation was examined as an example of a nonlinear 
problem treated by asymptotic methods. The analysis using linearly 
independent solutions at infinity allowed accurate computations of the 
solutions in the boundary layer and suggested that shooting methods applied 
to the problem should integrate backward from infinity rather than starting 
the solution at the origin. 
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Table 1. Asymptotic Approximations for Jg(x), Bessel Function 
of First Kind and Order 6 for Large x. 


X 

Exact 

value 

Newton's Method 

Approximation 
First Zeroth 

Eq. (2.33a) Eq. (2.34a) 

Perturbation Method 

Approximation 
Zeroth First 

Eq. (2.35a) Eq. (2.36a) 

5.98 

.243406659478 

0.00000000 

0.00000000 

-.15106856 

-1.01572984 

6.00 

.245836863364 

.77225975 

.22293219 

- .15679786 

-1.00738895 

6.20 

.268597150751 

.44343614 

.22807271 

-.20697376 

-.91224979 

6.40 

.289958360953 

.39840139 

.23792512 

- .24748793 

- .79352271 

6.60 

.309275672516 

.38345790 

.24952179 

- .27709881 

- .65697857 

6.80 

.325899479373 

.37851419 

.26122912 

-.29500488 

-.50844102 

7.00 

.339196604983 

.37695478 

.27182801 

- .30086304 

-.35363680 

7.20 

.348572602373 

.37556888 

.28032534 

- .29479002 

-.19805436 

7.40 

.353494359028 

.37245303 

.28589146 

- .27734706 

- .04681405 

7.60 

.353512167554 

.36635530 

.28783970 

- .24950883 

.09544723 

7.80 

.348280396189 

.35642320 

.28562343 

-.21261767 

.22467440 

8.00 

.337575900114 

.34209466 

.27884054 

-.16832536 

.33747340 

8.20 

.321313356102 

.32304740 

.26724063 

-.11852458 

.43116997 

8.40 

.299556779154 

.29917343 

.25073178 

- .06527280 

.50384658 

8.60 

.272526588542 

.27056419 

.22938517 

- .01071154 

.55435761 

8.80 

.240601729158 

.23749901 

.20343630 

.04301599 

.58232254 

9.00 

.204316517680 

.20043323 

.17328187 

.09384575 

.58809821 

9.20 

.164352066323 

.15998377 

.13947200 

.13986925 

.57273137 

9.40 

.121522333286 

.11691131 

.10269746 

.17940273 

.53789346 

9.60 

.076755051090 

.07209854 

.06377198 

.21104660 

.48579985 

9.80 

.031067984231 

.02652454 

.02361001 

.23373316 

.41911637 

10.00 

-.014458842085 

-.01876427 

-.01679964 

.24676089 

.34085607 

10.20 

-.058714127480 

- .06268601 

- .05642347 

.24981441 

.25426955 

10.40 

-.100588533517 

-.10415719 

- .09421547 

.24296933 

.16273215 

10.60 

-.139009092875 

-.14212759 

-.12915072 

.22668229 

.06963152 

10.80 

-.172973630937 

-.17561498 

- .16025935 

.20176666 

-.02174126 

11.00 

-.201584000874 

-.20373856 

-.18665967 

.16935504 

-.10829367 

11.20 

-.224076937511 

-.22574999 

-.20758919 

.13085024 

- .18722411 

11.40 

-.239851384024 

-.24106084 

-.22243262 

.08786670 

- .25610192 

11.60 

-.248491238197 

-.24926544 

-.23074553 

.04216483 

-.31293313 

11.80 

-.249782599151 

-.25015826 

-.23227284 

-.00441932 

- .35621012 

12.00 

-.243724767229 

- .24374494 

- .22696138 

- .05004610 

-.38494417 

12.20 

-.230534453434 

-.23024670 

-.21496589 

-.09294334 

- .39868021 

12.40 

-.210642883654 

-.21009751 

-.19664810 

-.13147387 

-.39749390 

12.60 

-.184685728838 

-.18393416 

-.17256878 

-.16419667 

-.38197150 

12.80 

-.153486046509 

-.15257939 

-.14347284 

-.18991972 

-.35317358 

13.00 

-.118030672130 

-.11701838 

-.11026792 

- .20774240 

- .31258420 

13.20 

-.079440741520 

-.07836949 

- .07399710 

- .21708612 

-.26204745 

13.40 

-.038937248376 

-.03785000 

-.03580641 

- .21771198 

- .20369389 

13.60 

.002197264605 

.00326205 

.00309151 

- .20972517 

- .13985934 

13.80 

.042659624317 

.04366858 

.04145673 

- .19356595 

- .07299893 

14.00 

.081168183426 

.08209341 

.07806274 

- .16998779 

-.00559950 

14.20 

.116504774373 

.11732408 

.11173727 

- .14002372 

.05990705 

14.40 

.147555004998 

.14825196 

.14140177 

- .10494224 

.12122613 

.14.60 

.173345475655 

.17390930 

.16610813 

- .06619474 

.17627638 
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Table 1. Continued. 


Newton's Method Perturbation Method 


Approximation Approximation 


X 

Exact 

value 

First 
Eq. (2.33a) 

Zeroth 
Eq. (2.34a) 

Zeroth 
Eq. (2.35a) 

First 
Eq. (2.36a) 

14.80 

.193076585993 

.19350189 

.18507128 

-.02535643 

.22325627 

15.00 

.206149737480 

.20643617 

.19769657 

.01593665 

.26069971 

15.20 

.212187920995 

.21233968 

.20360074 

.05605233 

.28751874 

15.40 

.211048901297 

.21107417 

.20262591 

.09342595 

.30303224 

15.60 

.202830464206 

.20274078 

.19484573 

.12662057 

.30697995 

15.80 

.187867468976 

.18767707 

.18056363 

.15438157 

.29952153 

16.00 

.166720737703 

.16644584 

.16030311 

.17568339 

.28122115 

17.00 

.000715333443 

.00028441 

.00027517 

.16920613 

.07047137 

18.00 

-.155956234195 

- .15618751 

-.15168840 

.01205744 

-.17431553 

19.00 

-.178767171544 

- .17870512 

- .17410631 

-.14737116 

- .24951527 

20.00 

-.055086049564 

- .05487232 

-.05360289 

-.16665635 

-.10973177 
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TABLE 2. 


r»4-1 n n 

Solution of Equation (3.19) c - [i (-dc)]/[d i (-dc)] where i (z) 

2 

is the nth Repeated Integral of the Error Function and d = 1/2. 


n 

P 

C 

5 

2.5 

.427 

4 

2.0 

. 473 

3 

1.5 

.537 

2 

1.0 

.639 

1 

0.5 

.840 

0 

0.0 

00 
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Table 3. Summary of Numerical Results for the Falkner-Skan Equation 
from the modified Newton iteration 


p 

c 

w (°°) 

f' ' (0) 

-w' ' (0) 

f' ' (0) 
Published 

Ref. 

No. 

Number of 
Iterations 

w o 

2.4 

.495 

- .4614 

1.8384 

1.837 

12 

7 

w o 

« - c 

2.4 

.495 

- .45985 

1.8378 

1.837 

12 

1 

w o 

- W< 7 > 

2.0 

.473 

- .497 

1.6873 
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